Estimating diarrhea cases prevented by WASH among children < 3y, indexed by relative wealth rank, per grid-cell and district
Pearl Ante-Testard
2023-03-09
Here, diarrhea cases prevented by combined WASH per wealth rank were calculated by:
Getting the absolute difference between diarrhea cases among intervention group under 3y and diarrhea cases among controls under 3y, indexed by wealth rank, where:
diarrhea cases among controls = children under 3y multiplied by the proportion of diarrhea among controls by wealth rank (which was estimated from the GAM) per grid-cell and aggregated by district
diarrhea cases among intervention = children under 3y multiplied by the proportion of diarrhea among intervention by wealth rank (estimated from the GAM) per grid-cell and aggregated by district
GAM estimates were based on the WASH Benefits Bangladesh cluster randomized trial. We used these GAM estimates to extrapolate to the whole country of Bangladesh, while masking for urban areas since the WASHB trial were conducted in rural areas in Gazipur, Mymensingh and Tangail districts.
library(here)
library(cowplot)
library(sp)
library(scales)
library(ggpubr)
library(patchwork)
source(here::here("3-secondary-analysis/R", "0-config.R"))# world population in 2014 (from worldpop)
pop_worldpop_bgd <- readRDS(file = here::here("1-data",
"1-temp",
"pop-under3-grid-cell-bgd.rds"))
# wealth index layer in 2011 (from worldpop)
wealth_worldpop_bgd <- readRDS(file = here::here("1-data",
"1-temp",
"wealthindex-grid-cell-bgd.rds"))
# diarrhea proportion among the controls estimated from GAM (WASHB trial)
diar_prev_control_wealth <- readRDS(file = here::here("1-data","1-temp",
"gam-diar-prevalence-control-wealthrank.RDS"))
# diarrhea proportion among the intervention group estimated from GAM (WASHB trial)
diar_prev_interv_wealth <- readRDS(file = here::here("1-data","1-temp",
"gam-diar-prevalence-interv-wealthrank.RDS"))
# shapefiles (from GADM)
BGD_shp_adm0 <- st_read(here::here("1-data", "0-gadm",
"gadm41_BGD_0.json"))## Reading layer `gadm41_BGD_0' from data source
## `/Users/pearlante/Library/CloudStorage/Box-Box/washb-sep-season-proj/washb-sep-season/washb-sep-season/1-data/0-gadm/gadm41_BGD_0.json'
## using driver `GeoJSON'
## Simple feature collection with 1 feature and 2 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 88.0106 ymin: 20.7411 xmax: 92.6737 ymax: 26.6341
## Geodetic CRS: WGS 84BGD_shp_adm2 <- st_read(here::here("1-data", "0-gadm",
"gadm41_BGD_2.json")) ## Reading layer `gadm41_BGD_2' from data source
## `/Users/pearlante/Library/CloudStorage/Box-Box/washb-sep-season-proj/washb-sep-season/washb-sep-season/1-data/0-gadm/gadm41_BGD_2.json'
## using driver `GeoJSON'
## Simple feature collection with 64 features and 13 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 88.0113 ymin: 20.7411 xmax: 92.6737 ymax: 26.6341
## Geodetic CRS: WGS 84# urban raster layer (from Francois based on the Global Human Settlement)
ras <- readRDS(here::here("1-data", "Urban_BGD"))
# WASH Benefits Bangladesh clusters
df_gps <- read_dta(file = here::here("1-data", "0-untouched",
"6. WASHB_Baseline_gps.dta"))
df_publicid <- read_csv(file = here::here("1-data", "0-untouched",
"public-ids.csv"))
# individual datasets included in the analysis
df_analysis_svy12 <- readRDS(file = here::here("1-data", "2-final",
"enrol_diar_tr_wealth_indiv_svy012.rds")) %>%
filter(svy!=0)# rounding wealth rank to 1 decimal point before merging to other datasets
wealth_worldpop_bgd <- wealth_worldpop_bgd %>%
mutate(wealth_rank=round(wealth_rank, 1))
diar_prev_control_wealth <- diar_prev_control_wealth %>%
mutate(wealth_rank=round(wealth_rank, 1))
diar_prev_interv_wealth <- diar_prev_interv_wealth %>%
mutate(wealth_rank=round(wealth_rank, 1)) ## join worldpop population 2014 and worldpop wealth 2011
wealth_popunder3_gps <- st_join(pop_worldpop_bgd, wealth_worldpop_bgd) %>%
filter(!is.na(bgd2011wipov.tif))
wealth_popunder3_gps <-wealth_popunder3_gps[!duplicated(wealth_popunder3_gps$bgd_ppp_2014_1km_Aggregated.tif), ]
## join with GAM proportion estimates for the control group
diar_wealth_pop_gps_control <- inner_join(wealth_popunder3_gps %>% as.data.frame(),
diar_prev_control_wealth,
by="wealth_rank")
diar_wealth_pop_gps_control_sf <- st_as_sf(diar_wealth_pop_gps_control)
diar_wealth_pop_gps_control_sf <-diar_wealth_pop_gps_control_sf[
!duplicated(diar_wealth_pop_gps_control_sf$bgd_ppp_2014_1km_Aggregated.tif), ]
diar_wealth_pop_gps_control_sf <- diar_wealth_pop_gps_control_sf %>%
mutate(pop_under3_diar_wealth_control=pop_3andunder*diar7d)
## join with GAM proportion estimates for the intervention group
diar_wealth_pop_gps_interv <- inner_join(wealth_popunder3_gps %>% as.data.frame(),
diar_prev_interv_wealth,
by="wealth_rank")
diar_wealth_pop_gps_interv_sf <- st_as_sf(diar_wealth_pop_gps_interv)
diar_wealth_pop_gps_interv_sf <-diar_wealth_pop_gps_interv_sf[
!duplicated(diar_wealth_pop_gps_interv_sf$bgd_ppp_2014_1km_Aggregated.tif), ]
diar_wealth_pop_gps_interv_sf <- diar_wealth_pop_gps_interv_sf %>%
mutate(pop_under3_diar_wealth_interv=pop_3andunder*diar7d)
diar_wealth_pop_gps_interv_sf <- diar_wealth_pop_gps_interv_sf[,(c("bgd_ppp_2014_1km_Aggregated.tif","bgd2011wipov.tif",
"wealth_rank","Arms","geometry",
"pop_under3_diar_wealth_interv"))]
## merge control and intervention df
diar_wealth_pop_gps_interv_control_sf <- inner_join(diar_wealth_pop_gps_control_sf,
diar_wealth_pop_gps_interv_sf %>% as.data.frame(),
by=c("bgd_ppp_2014_1km_Aggregated.tif",
"bgd2011wipov.tif",
"wealth_rank"))
## create the excess diarrhea between the intervention and control groups
diar_wealth_pop_gps_interv_control_sf <- diar_wealth_pop_gps_interv_control_sf %>%
mutate(excess_diar= abs(pop_under3_diar_wealth_interv-pop_under3_diar_wealth_control))
## join diar_wealth_pop_gps_interv_control_sf with the shapefile (needed for district estimates)
#bgd_shp_adm2 <- st_make_valid(BGD_shp_adm2)
#diar_wealth_pop_gps_interv_control_sf <- st_make_valid(diar_wealth_pop_gps_interv_control_sf)
#joined <- st_join(bgd_shp_adm2, diar_wealth_pop_gps_interv_control_sf)
## aggregate excess_diar by district
#joined_grouped <- joined %>%
#group_by(NAME_2) %>%
#mutate(excess_diar_sum = sum(excess_diar, na.rm = T),
# excess_diar_per1000under3 = excess_diar_sum/1000)## transforming urban/rural layer to dataframe
r.spdf <- as(ras, "SpatialPixelsDataFrame")
r.df <- as.data.frame(r.spdf)
head(r.df)## layer x y
## 1 0 88.40467 26.62632
## 2 0 88.40467 26.61796
## 3 0 88.41517 26.61796
## 4 0 88.39417 26.60960
## 5 0 88.40467 26.60960
## 6 0 88.39417 26.60124# urban (needed for the urban map)
r.df.urban <- subset(r.df, layer==1)
## transforming into points
r.df.sf <- r.df %>%
st_as_sf(coords = c("x","y"), crs = 4326)
## merging the coordinates into the df
joined.r <- st_join(r.df.sf, diar_wealth_pop_gps_interv_control_sf)
## making the urban layers into NAs for the maps in grid-cell
joined.r <- joined.r %>%
filter(!is.na(excess_diar))
joined.r$excess_diar[joined.r$layer == 1] <- NA
## join the shapefile with the layers for the district estimates
bgd_shp_adm2 <- st_make_valid(BGD_shp_adm2)
diar_wealth_pop_gps_interv_control_sf <- st_make_valid(diar_wealth_pop_gps_interv_control_sf)
joined <- st_join(bgd_shp_adm2, joined.r)
#joined_grouped_urban <- subset(joined, layer==1)
#joined_grouped_rural <- subset(joined, layer==0)
#joined <- joined %>%
# filter(!is.na(excess_diar))
#joined_grouped_rural <- joined %>%
# group_by(NAME_2) %>%
#mutate(excess_diar_sum = sum(excess_diar, na.rm = T),
#excess_diar_per1000under3 = excess_diar_sum/1000)df_gps_publicid <- inner_join(df_gps, df_publicid, by = "dataid")
df_gps_publicid <- df_gps_publicid[,-c(1,5,6)]
colnames(df_gps_publicid) <- c("qgpslong","qgpslat","entrydate","clusterid",
"dataid","block")
df_gps_publicid_washb_sepseason <- inner_join(df_gps_publicid, df_analysis_svy12, by = c("dataid",
"clusterid",
"block"))
df_gps_publicid_clust <- df_gps_publicid_washb_sepseason %>%
st_as_sf(coords = c("qgpslong", "qgpslat"), crs = 4326) %>%
group_by(clusterid) %>%
dplyr::summarise() %>%
st_cast("POLYGON")
df_gps_publicid_block <- df_gps_publicid_washb_sepseason %>%
st_as_sf(coords = c("qgpslong", "qgpslat"), crs = 4326) %>%
group_by(block) %>%
dplyr::summarise() %>%
st_cast("POLYGON")WASH Benefits Bangladesh clusters
map_bgd_wash_clust <- ggplot(df_gps_publicid_clust) +
geom_sf(data = BGD_shp_adm0, fill = "grey") +
#geom_sf(color="blue") +
geom_point(size=0.5, shape = 21, aes(geometry=geometry),
stat = "sf_coordinates") +
geom_sf(data = BGD_shp_adm2, fill = NA) +
theme_void() +
#theme(legend.position = "bottom") +
ggtitle("WASH Benefits Bangladesh trial clusters included in the analysis") +
theme(plot.title = element_text(hjust = 0.5))
map_bgd_wash_clust## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data
# the trial was located in Gazipur, Mymensingh and Tangail districts WASH Benefits Bangladesh block
map_bgd_wash_block <- ggplot(df_gps_publicid_block) +
geom_sf(data = BGD_shp_adm0, fill = "grey") +
#geom_sf(color="blue") +
geom_point(size=1, shape = 21, aes(geometry=geometry),
stat = "sf_coordinates") +
geom_sf(data = BGD_shp_adm2, fill = NA) +
theme_void() +
#theme(legend.position = "bottom") +
ggtitle("WASH Benefits Bangladesh trial blocks included in the analysis") +
theme(plot.title = element_text(hjust = 0.5))
map_bgd_wash_block## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data
# the trial was located in Gazipur, Mymensingh and Tangail districts Population density and WASH Benefits Study clusters
# Colors
pal <- hcl.colors(10, "mako", rev = TRUE)
# Plot Density
bgd_density <- ggplot() +
geom_sf(data = pop_worldpop_bgd,
aes(fill = popdensity),
color=NA,
show.legend=TRUE) +
scale_fill_gradientn(colors=pal,
values=rescale(seq(0,10,by=2)),
limits=c(0,10),
#n.breaks = 8
breaks = c(0,2,4,6,8,10)) +
geom_point(data=df_gps_publicid_clust, size=0, shape = 21, color = "#E34D34" ,
aes(geometry=geometry), stat = "sf_coordinates") +
geom_sf(data = bgd_shp_adm2, fill = NA, color="black") +
coord_sf() +
theme_void() +
theme(legend.position = "right",text = element_text(size=7.5)) +
labs(
fill = "Population density of \nchildren < 3y per grid-cell")
bgd_density <- bgd_density + annotate(
geom = "curve", x=91, y=25.8, xend = 90.5, yend = 25,
curvature = .1, arrow = arrow(length = unit(2, "mm")), col="#E34D34"
) +
annotate("text", x=91.2, y=25.9, label= "WASH Benefits Study clusters (n=360)",
col="#E34D34", size=2.5, parse=F)
bgd_density## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data
Wealth index layer
# Read tif file for wealth index, 2011
wipov_tif=read_stars(here::here("1-data", "0-worldpop",
"bgd2011wipov.tif"))
wipov_sf=st_as_sf(wipov_tif)
# Wealth rank
wipov_sf <- wipov_sf %>%
mutate(wealth_rank = rank(bgd2011wipov.tif)/nrow(.),
wealth_rank = round(wealth_rank,1))
# Colors
pal <- hcl.colors(10, "Rocket", rev = T)
bgd_wealth_cont <- ggplot() +
geom_sf(data = wipov_sf,
aes(fill = wealth_rank),
color=NA,
show.legend=TRUE) +
scale_fill_gradientn(colors=pal,
values=rescale(seq(0,1,by=0.2)),
limits=c(0,1),
#n.breaks = 8
breaks = seq(0,1,by=0.2)) +
geom_sf(data=BGD_shp_adm2, fill=NA, color="black") +
coord_sf() +
theme_void()+
theme(legend.position = "right",text = element_text(size=7.5)) +
labs(
fill = "Relative rank based \non the wealth index")
bgd_wealth_cont
#ggsave(bgd_wealth_cont, file = here::here("3-secondary-analysis",
# "output",
# "map-wealthrank-bgd.png"), height=8,width=7)Masked urban layer
The urban/rural layer from Francois Rerolle was based on the Global Human Settlement.
# urban/rural layer from Francois
#plot_ras <- plot(ras)
plot_maskedurban <- ggplot() +
geom_sf(data=BGD_shp_adm2, color="black", fill="grey") +
geom_tile(data=r.df.urban, aes(x=x,y=y), fill="darkblue") +
coord_sf() +
theme_void() +
labs(subtitle = "Urban areas that were masked in \nthe spatial analysis (blue color)") +
#ggtitle("Urban areas (blue color) that were masked in the spatial analysis") +
theme(text = element_text(size=7.5), plot.caption = element_text(hjust = 0.5))
plot_maskedurban
Diarrhea cases prevented for children <3y per grid-cell at a 1 km res.
#joined.r$excess_diar[joined.r$layer == 1] <- NA
#Min. 1st Qu. Median Mean 3rd Qu. Max.
# 0.00 1.18 2.07 2.69 3.58 27.31
colnames(joined.r) <- c("layer","bgd_ppp_2014_1km_Aggregated.tif","pop_0_14","pop_3andunder","popdensity",
"bgd2011wipov.tif","wealth_rank","wealthdensity",".idx","Arms_control","block",
"dummy","diar7d","SE","CI_upper","CI_lower","pop_under3_diar_wealth_control",
"Arms_interv","pop_under3_diar_wealth_interv","excess_diar","geometry_point","geometry")
# layer
# 0 1
#111519 26870
joined.r_2 <- joined.r[,-c(21)]
## plot the diarrhea cases prevented by WASH
pal <- hcl.colors(15, "Plasma", rev = T)
plot_diar_averted_bgd_gridcell_rural <- ggplot() +
#geom_tile(data=joined_rural_2, aes(x=x,y=y), fill=NA) +
geom_sf(data = joined.r,
aes(color=excess_diar*4),
#color=NA,
show.legend=TRUE) +
scale_color_gradientn(colors=pal,
values=rescale(seq(0,28,by=4)),
limits=c(0,28),
n.breaks = 8,
breaks = seq(0,28,by=4),
na.value = 'darkgrey'
) +
#scale_y_continuous(breaks=seq(0,28,by=4)) +
#geom_tile(data=joined_urban.r, aes(x=x,y=y), fill="grey") +
geom_sf(data=BGD_shp_adm2, fill=NA, color="black") +
coord_sf() +
#new_scale_fill() +
theme_void()+
theme(legend.position = "right",text = element_text(size=7.5)) +
labs(
color = "Diarrhea cases prevented for \nchildren < 3y per month \nper grid-cell")
plot_diar_averted_bgd_gridcell_rural <- plot_diar_averted_bgd_gridcell_rural +
annotate(
geom = "curve", x=91, y=21.5, xend = 91.8, yend = 22.3,
curvature = .1, arrow = arrow(length = unit(1.75, "mm")), col="darkgrey"
) +
annotate("text", x=91, y=21.3, label= "Masked urban areas \n(excluded from estimates)",
col="darkgrey", size=2.5, parse=F)
plot_diar_averted_bgd_gridcell_rural
#ggsave(plot_diar_averted_bgd_gridcell_rural, file = here::here("3-secondary-analysis",
# "output",
# "map-diar-under3-wealthrank-bgd-rural-gridcell-2.jpg"), height=8,width=7)
#0,4,8,12,16,20,24,28,32Diarrhea cases prevented per 1000 children < 3y per week in each district
## subset rural layers by district
joined_grouped_rural <- subset(joined, layer==0)
joined_grouped_rural <- joined_grouped_rural %>%
group_by(NAME_2) %>%
mutate(excess_diar_sum = sum(excess_diar, na.rm = T),
excess_diar_per1000under3 = excess_diar_sum/1000)
## plot the diarrhea cases prevented by WASH excluding urban areas
#pal <- hcl.colors(5, "Viridis", rev = F, alpha = 0.9)
#my_breaks = c(0, 1, 2, 3, 4)
joined_grouped_rural <-joined_grouped_rural[
!duplicated(joined_grouped_rural$NAME_2), ]
plot_diar_averted_bgd_rural_dist <- ggplot() +
geom_sf(data = joined_grouped_rural,
aes(fill=excess_diar_per1000under3*4),
color=NA,
show.legend=TRUE) +
#scale_fill_gradientn(colors = pal, guide="legend") +
scale_fill_viridis_c(alpha = 1,begin = 0,end = 1,aesthetics = "fill",direction = -1,
values=scales::rescale(seq(0,20,by=2)),
breaks = seq(0,18,by=2),limits=c(0,20)
)+
#scale_y_continuous(breaks=seq(0,20,by=4)) +
geom_sf(data=BGD_shp_adm2, fill=NA, color="black") +
coord_sf() +
#new_scale_fill() +
theme_void()+
theme(legend.position = "none",text = element_text(size=7.5)) #+
#labs(
#fill = "Diarrhea cases prevented \nper 1000 children < 3y per \nmonth in each district \n(masked urban areas)")
plot_diar_averted_bgd_rural_dist_print <- ggplot() +
geom_sf(data = joined_grouped_rural,
aes(fill=excess_diar_per1000under3*4),
color=NA,
show.legend=TRUE) +
#scale_fill_gradientn(colors = pal, guide="legend") +
scale_fill_viridis_c(alpha = 1,begin = 0,end = 1,aesthetics = "fill",direction = -1,
values=scales::rescale(seq(0,20,by=2)),
breaks = seq(0,18,by=2),
limits=c(0,20)
)+
#scale_y_continuous(breaks=seq(0,18,by=2)) +
geom_sf(data=BGD_shp_adm2, fill=NA, color="black") +
coord_sf() +
#new_scale_fill() +
theme_void()+
theme(legend.position = "right",text = element_text(size=7.5)) +
labs(
fill = "Diarrhea cases prevented \nper 1000 children < 3y per \nmonth in each district")
plot_diar_averted_bgd_rural_dist_print 
#ggsave(plot_diar_averted_bgd_rural_dist, file = here::here("3-secondary-analysis",
# "output",
# "map-diarprevented-maskedurban-bgd-2.png"), height=8,width=7)Barplot
## barplot of the diarrhea cases prevented per district
#joined_grouped_rural <-joined_grouped_rural[
#!duplicated(joined_grouped_rural$NAME_2), ]
joined_grouped_rural <- joined_grouped_rural %>%
mutate(NAME_2 = factor(NAME_2,levels=joined_grouped_rural$NAME_2[order(joined_grouped_rural$excess_diar_per1000under3)]))
options(scipen=999)
#pal <- hcl.colors(5, "Viridis", rev = F, alpha = 0.9)
barplot_diar_averted_bgd_rural_dist <- ggplot(joined_grouped_rural, aes(x=NAME_2, y=excess_diar_per1000under3*4, fill=excess_diar_per1000under3*4)) +
geom_bar(stat = "identity") +
scale_fill_viridis_c(alpha = 1,begin = 0,end = 1,aesthetics = "fill",direction = -1,
values=scales::rescale(seq(0,20,by=2)),
#breaks = seq(0,18,by=2),
limits=c(0,20)
)+
scale_y_continuous(breaks=seq(0,18,by=2)) +
geom_hline(yintercept=seq(0,20,by=2), colour="white") +
coord_flip() +
xlab(" ") +
#scale_y_continuous(position = 'right') +
ylab(
"Diarrhea cases prevented per 1000 children < 3y per month \nin each district") +
theme_bw() +
theme(text = element_text(size=7.5)) +
theme(legend.position = "none") +
theme(panel.grid.major=element_blank(), panel.grid.minor=element_blank())
barplot_diar_averted_bgd_rural_dist
#ggsave(barplot_diar_averted_bgd_rural_dist, file = here::here("3-secondary-analysis",
# "output",
# "barplot-diarprevented-maskedurban-bgd-2.png"), height=8,width=7)Combined plots
map_district_composite <- barplot_diar_averted_bgd_rural_dist + inset_element(plot_diar_averted_bgd_rural_dist, left = 0.4, bottom = 0, right = 1, top = 0.9)
map_district_composite
plot_composite <- ggarrange(bgd_density, bgd_wealth_cont, plot_diar_averted_bgd_gridcell_rural,
map_district_composite, #+ rremove("x.text"),
labels = c("A", "B", "C", "D"),
ncol = 2, nrow = 2, align = "hv", widths = c(2, 2))## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data## Warning: Graphs cannot be vertically aligned unless the axis parameter is set.
## Placing graphs unaligned.## Warning: Graphs cannot be horizontally aligned unless the axis parameter is set.
## Placing graphs unaligned.plot_composite
ggsave(plot_composite, file = here::here("3-secondary-analysis",
"output",
"map-composite-bgd-17.png"), height=10,width=10)
#plot_composite_2 <- plot_grid(plot_maskedurban, plot_diar_averted_bgd_gridcell_rural, align = "hv",
# labels = c("A", "B"))
#plot_composite_2
#ggsave(plot_composite_2, file = here::here("3-secondary-analysis",
# "output",
# "map-composite-masked-bgd-2.png"), height=8,width=8)# plot
#ggsave(plot_composite, file = here::here("3-secondary-analysis",
# "output",
# "map-diar-under3-wealthrank-bgd-rural.png"), height=8,width=7)
#ggsave(plot_composite_2, file = here::here("3-secondary-analysis",
# "output",
# "map-diar-under3-wealth-gridcell-bgd-rural.png"), height=8,width=7)
# data frame
#saveRDS(joined_grouped, file = here::here("1-data",
# "1-temp",
# "diar-wealthrank-cases-prevented-district-v3.rds"))
#write_csv(joined_grouped, file = here::here("1-data",
# "1-temp",
# "diar-wealthrank-cases-prevented-district-v3.csv"))sessionInfo()## R version 4.2.1 (2022-06-23)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur ... 10.16
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] raster_3.6-11 colorspace_2.0-3 ggrepel_0.9.1 stars_0.5-6
## [5] abind_1.4-5 viridis_0.6.2 viridisLite_0.4.1 forcats_0.5.2
## [9] stringr_1.4.1 dplyr_1.0.10 purrr_0.3.5 readr_2.1.2
## [13] tidyr_1.2.1 tibble_3.1.8 tidyverse_1.3.2 sf_1.0-8
## [17] tidygam_0.1.0.9000 tidymv_3.3.2 mgcv_1.8-41 nlme_3.1-157
## [21] haven_2.5.1 patchwork_1.1.2 ggpubr_0.4.0 ggplot2_3.4.0
## [25] scales_1.2.1 sp_1.5-0 cowplot_1.1.1 here_1.0.1
##
## loaded via a namespace (and not attached):
## [1] googledrive_2.0.0 ggsignif_0.6.3 ellipsis_0.3.2
## [4] class_7.3-20 rprojroot_2.0.3 fs_1.5.2
## [7] rstudioapi_0.14 proxy_0.4-27 farver_2.1.1
## [10] bit64_4.0.5 fansi_1.0.3 lubridate_1.8.0
## [13] xml2_1.3.3 codetools_0.2-18 splines_4.2.1
## [16] cachem_1.0.6 knitr_1.40 jsonlite_1.8.3
## [19] broom_1.0.1 dbplyr_2.2.1 compiler_4.2.1
## [22] httr_1.4.4 backports_1.4.1 assertthat_0.2.1
## [25] Matrix_1.5-1 fastmap_1.1.0 gargle_1.2.1
## [28] cli_3.4.1 s2_1.1.0 htmltools_0.5.3
## [31] tools_4.2.1 gtable_0.3.1 glue_1.6.2
## [34] wk_0.6.0 Rcpp_1.0.9 carData_3.0-5
## [37] cellranger_1.1.0 jquerylib_0.1.4 vctrs_0.5.0
## [40] lwgeom_0.2-8 xfun_0.33 rvest_1.0.3
## [43] lifecycle_1.0.3 rstatix_0.7.0 googlesheets4_1.0.1
## [46] terra_1.6-47 vroom_1.5.7 ragg_1.2.3
## [49] hms_1.1.2 parallel_4.2.1 yaml_2.3.6
## [52] gridExtra_2.3 sass_0.4.2 stringi_1.7.8
## [55] highr_0.9 e1071_1.7-11 systemfonts_1.0.4
## [58] rlang_1.0.6 pkgconfig_2.0.3 evaluate_0.17
## [61] lattice_0.20-45 labeling_0.4.2 bit_4.0.4
## [64] tidyselect_1.2.0 magrittr_2.0.3 R6_2.5.1
## [67] generics_0.1.3 DBI_1.1.3 pillar_1.8.1
## [70] withr_2.5.0 units_0.8-0 modelr_0.1.9
## [73] crayon_1.5.2 car_3.1-0 KernSmooth_2.23-20
## [76] utf8_1.2.2 tzdb_0.3.0 rmarkdown_2.16
## [79] grid_4.2.1 readxl_1.4.1 reprex_2.0.2
## [82] digest_0.6.30 classInt_0.4-7 textshaping_0.3.6
## [85] munsell_0.5.0 bslib_0.4.0